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Summary 


During the current reporting period, January 1 to December 31, 1989, general problems 
associated with on-board trajectory optimization, propulsion system cycle selection, and with the 
synthesis of guidance laws were addressed for ascent to low-Earth-orbit of an air-breathing, 
single-stage-to-orbit vehicle. This report follows a previous one entitled "Trajectory Optimization 
and Guidance Law Development for National Aerospace Plane Applications" and dated December 
1988. The work reported herein builds directly upon the analytical results presented in that 
document. A good portion of this work focused on making improvements to the vehicle models 
employed. The NASA "Generic Hypersonic Aerodynamic Model Example" and the "Langley 
Accelerator" aerodynamic data sets were acquired and implemented. Work pertaining to the 
development of purely analytic aerodynamic models also continued at a low level. A generic model 
of a multi-mode propulsion system was developed that includes turbojet, ramjet, scramjet, and 
rocket engine cycles. Provisions were made in the dynamic model for a component of thrust 
normal to the flight path. Computational results, which characterize the nonlinear sensitivity of 
scramjet performance to changes in vehicle angle of attack, were obtained and incorporated into the 
engine model. Additional trajectory constraints were also introduced. The constraints now treated 
are: maximum dynamic pressure, maximum aerodynamic heating rate per unit area, angle of attack 
and lift limits, and limits on acceleration both along and normal to the flight path. 

The remainder of the research effort focused, for the most part, on required modifications to the 
previously derived algorithm when the model complexity cited above was added. In particular, 
analytic switching conditions were derived which, under appropriate assumptions, govern optimal 
transition from one propulsion mode to another for two cases: the case in which engine cycle 
operations can overlap, and the case in which engine cycle operations are mutually exclusive. The 
resulting guidance algorithm was implemented in software and exercised extensively. It was found 
that the approximations associated with the assumed time scale separation employed in this work 
are reasonable except over the Mach range from roughly 5 to 8. This phenomenon is due to the 
very large thrust capability of scramjets in this Mach regime when sized to meet the requirement for 
ascent to orbit. Very little mass penalty is induced by the resulting inaccuracies in the trajectory 
over this region because it is traversed rapidly. However, the reduced solution climb paths prove to 
be unfeasible within this Mach range when subject to the full model dynamics and active trajectory 
constraints. These difficulties were successfully overcome by accounting for flight path angle and 
flight path angle rate in construction of the flight path over this Mach range. The resulting 
algorithm provides the means for rapid near-optimal trajectory generation and propulsion cycle 
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selection over the entire Mach range from take-off to orbit given a realistic nonlinear vehicle model 
and all pertinent trajectory constraints. 

The only significant problem area encountered to date relates to the lack of a general theory for 
singularly perturbed systems that are. subject to state-variable inequality constraints. Such 
constraints are common to a wide class of flight vehicles but have received little attention in the 
literature when the dynamic system is singularly perturbed. A study was initiated in this area and it 
was found that, when the reduced solution lies on a state-variable inequality constraint boundary, 
the boundary layer trajectories are of finite time in the stretched time scale. The possibility of 
costate discontinuites at the juncture between constrained and unconstrained arcs makes direct 
application of existing theory difficult at best. A transformation technique was identified that 
eliminates some of these difficulties, but at the cost of possibly increased system order and the 
introduction of singular arcs. Much work remains to be done in this area. 

Work on development of simple, efficient algorithms for prediction of vehicle aerodynamic and 
propulsive performance have continued during the present phase of the program. Improvements in 
modeling of the hypersonic lifting body module have eliminated previous discrepancies between 
measured and predicted aerodynamic behavior. Several modes of data entry are now implemented 
making assessment of a given vehicle configuration very simple. An interactive program mode has 
been devised that makes possible direct and immediate assessment of configuration changes on 
selected vehicle performance paramaters. The algorithms developed in this program are of potential 
use in applications beyond those originally envisioned. 

Four conference papers have now been published which discuss most of the results of this 
research effort. A Ph.D. Dissertation that details the entire effort to date was published in 
Decemeber of 1989. A full-length paper entitled "Rapid Near-Optimal Trajectory Generation for 
Single-Stage-to-Orbit Airbreathing Vehicles" has been submitted for publication in the AIAA 
Journal of Guidance, Control and Dynamics and a new paper is now being prepared for the 1990 
AIAA GN&C Conference on the issue of state contraints in singularly perturbed systems. 


iii 


Table of Contents 


Page 

List of Figures v 

List of Symbols vi 

Section 1 - Introduction . 1 

Section 2 - General Problem Formulation 3 

Section 3 - Singular Perturbation Analysis 7 

3.1 Reduced Solution 7 

Lift as a Control 7 

Angle of Attack as a Control 10 

Bank Angle as a Control 12 

3.2 Boundary Layer Analysis 13 

Feedback Linearization - Lift as a Control 13 

Feedback Linearization - Angle of Attack as a Control 14 

Section 4 - Vehicle Models 16 

Section 5 - Numerical Results 17 

Section 6 - Conclusions and Recommendations 30 

6.1 Conclusions 30 

6.2 Recommendations 30 

6.3 Publications 31 

References 32 

Appendix A - State Inequality Constrained Boundary Layers 36 

Appendix B - Performance Modeling of Hypersonic Vehicles 52 


IV 



List of Figures 


Figure Eags 

1 . Partial Force and Moment Diagram 21 

2. Fuel specific impulse versus Mach number for the various propulsion cycle 

available to Vehicle Models 2-4 22 

3 . Reduced solution for Vehicle Model 1 as the maximum aerodynamic heating rate 

constraint (shown in Watts/cm2) is varied 23 

4 . Reduced solution for Vehicle Model 4 as the maximum axial acceleration limit 


(shown in g's) is varied 24 

5. Assumed variation in scramjet thrust with angle of attack 25 

6 . Reduced solution for Vehicle Model Number 3 - the effect of thrust variation with 

angle of attack is shown as the design angle of attack is varied 26 

7 . Reduced solution for Vehicle Model Number 2 - the effect of thrust variation with 

angle of attack is shown for a design angle of attack of 3 degrees 27 

8. Guided solution for Vehicle Model 4 - altitude time history 28 


9. Reduced solution for Vehicle Model 4 - comparison of the trajectory when flight 
path angle and flight path angle rate are included in the calculation of lift with the 
trajectory when they are not included , - 

A.l Illustration of finite costate rate at juncture between boundary layer and reduced 


solution trajectories 43 

B.l Hypersonic Airfoil 55 

B.2 Definition of Fuselage Reference Plane 56 

B.3 Definition of Body Shape Functions 58 



List of Symbols 

Cj = thrust specific fuel consumption for engine type j 
di = inequality constraint, dynamic pressure 
C2 = inequality constraint, aerodynamic heating 
C3 = inequality constraint, tangential acceleration 
C4 = inequality constraint, normal acceleration 
C Do = zero-lift drag coefficient 
CLa = lift curve slope 

D = aerodynamic drag 

E = mass specific energy 

Fc = thrust component aligned with the velocity vector 
Fs = thrust component normal to the velocity vector 
f = total fuel flow rate 

fmaxi = combined fuel flow rate for any number of independent engines of type i when operating 
at a stoichiometric fuel-to-air ratio 
g = acceleration due to gravity at mean sea level 

h = altitude, given by r - te 

H = Hamiltonian function 

J = performance index 

K = coefficient of the lift induced drag component 
L = aerodynamic lift 

m = vehicle mass 

M = Mach number 

n = total number of different engine cycles employed 

n-p = number of engine types employed for which fuel flow rate varies in direct proportion to (p 

but thrust varies in a nonlinear fashion with (p 
ni = component of acceleration aligned with the velocity vector 
n2 = component of acceleration normal to the velocity vector 

p = number of engine types employed that exhibit linear dependence of both fuel flow rate and 
thrust on q 

q = dynamic pressure 

Q = aerodynamic heating rate per unit area 

r = radial distance from the center of the Earth 
r£ = mean sea level 

s = aerodynamic reference area 

Sj = switching function for bang-bang control of qj 

t = time 

tf = final time 

T = net thrust 

Tk = combined thrust of any number of independent engines of type k, k=l to p 
U = pseudo control variable 

V = velocity 

W = Vehicle weight, taken at mean sea level (W = mg) 

a = angle of attack 

azL = angle of attack for zero lift 

y = flight path angle 

Se = elevon deflection 


VI 


e = perturbation parameter, 0 or 1 

exit = angle between T* and the body longitudinal axis 

T|j = throttle control, engine type j 

X x = costate, subscript denotes related state 

ji = gravitational constant for the Earth 

p = atmospheric density 

a = bank angle 

x = transformed time variable, x = t/e 

<Pi = fuel equivalence ratio, engine type i 

it = vector containing control for each engine type 


vii 


SECTION 1 
Introduction 

Emerging technology in many engineering fields, including hypersonic air-breathing 
propulsion, computational fluid dynamics, and high temperature materials, may soon make 
possible a vehicle configuration that has been the subject of study for over four decades 1 . This 
vehicle concept is commonly referred to as an aerospace plane. Its development, in one version or 
another, is being pursued by a number of industrialized nations. The current U.S. concept consists 
of a single-stage vehicle propelled, for the most part, by airbreathing engines. Most notable among 
the airbreathing cycles to be employed is that of the supersonic combustion ramjet or "scramjet. 
This aircraft is to be fueled by liquid hydrogen and will take-off and land horizontally on 
conventional runways. Operational objectives include hypersonic cruise in the upper atmosphere 
for long durations and the ability to accelerate to orbital velocity. Potential missions for such a 
vehicle include transportation to low-Earth-orbit, intercontinental passenger transportation, and a 
wide range of defense missions. This research effort is focused upon the particular mission of 
single-stage- to-orbit which promises, by the use of air-breathing hypersonic propulsion and 
greatly reduced launch operations, an order of magnitude reduction in the cost of placing payloads 
in low Earth orbit 2 - 3 . 

Even with the greatly improved fuel efficiency of airbreathing propulsion over current rocket 
engine technology, the ability to attain orbit in a single-stage vehicle will be marginal at best 4 . 
Trajectory optimization will play an important role in mission success for this reason. In fact, 
because the airbreathing propulsion system characteristics are sensitive to vehicle attitude and 
atmospheric conditions, precise trajectory control will be required. State-of-the-art launch vehicle 
guidance technology is heavily reliant on pre-mission, ground-based trajectory 
generation/optimization. In order to be cost effective, aerospace plane operations will have to 
approach those of modem commercial airlines. Technology dependent upon pre-mission, ground- 
based trajectory optimization is inadequate for this task; real-time, on-board trajectory optimization 
will be required 5 . 

The state of the art in trajectory optimization for complex nonlinear systems consists of a 
number of well developed numerical methods of solution. Unfortunately, these algorithms are 
poorly suited for on-board, real-time implementation. They are, in general, computationally 
intense, require an initial guess of the solution, and are lacking in global convergence 
characteristics. While some success in designing a reliable algorithm to numerically solve a two 
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point boundary value problem in an on-board computer has been achieved for orbit transfer 6 , the 
diverse mission requirements and complex control structure of a general purpose aerospace plane 
will likely require that structured methods for order reduction be employed. 

Energy state approximations and singular perturbation methods have been successfully 
employed to derive near-analytic trajectory optimization algorithms in the past. Near-optimal 
feedback guidance laws have also been obtained. These methods also contribute considerable 
insight into the nature of the optimal profiles and their relation to vehicle aerodynamic and 
propulsion characteristics. Early studies were devoted to fighter aircraft performance optimization 7 ' 
9. However, many of the modeling approximations employed for analysis of subsonic and 
supersonic aircraft optimal trajectories are not valid for a vehicle with hypersonic cruise and orbital 
capabilities. 

This research report presents an analysis of the problem of fuel-optimal ascent to low-Earth- 
orbitof an airbreathing, single-stage-to-orbit vehicle. Section II presents the problem formulation. 
A generic multi-mode propulsion system is defined which incorporates turbojet, ramjet, scramjet, 
and rocket engines. Inequality constraints on dynamic pressure, aerodynamic heating rate, and 
vehicle acceleration are also introduced. In Section HI an algorithm for generating fuel-optimal 
climb profiles is derived employing an energy state approximation. This algorithm results from 
application of the minimum principle to a low order dynamic model that includes general functional 
dependence on angle of attack and a normal component of thrust. Switching conditions are derived 
which, under appropriate assumptions, govern optimal transition from one propulsion mode to 
another. The use of bank angle to modulate the magnitude of the vertical component of lift is also 
investigated. A nonlinear transformation technique is employed to derive a feedback controller for 
tracking the computed trajectory. Section IV provides an overview of the vehicles models 
employed in this work. Section V provides a presentation and discussion of representative 
numerical results, and Section VI states conclusions drawn from this work. The main body of the 
report is followed by two appendices. Appendix A details an initial investigation into the 
characteristics of boundary layer systems when the reduced solution lies on a state-variable 
inequality constraint boundary. Appendix B details work performed in analytical vehicle model 
development 
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SECTION 2 

General Problem Formulation 

Consider atmospheric flight of a point mass over a spherical non-rotating Earth. The equations 
governing such flight can be reduced to a four state model as follows. 


* V(F C - D) 

E = — 


(1) 

m 


m = -f (r,E,7t,a) 


(2) 

• (F s + L) cosa 
67 ~ mV 

(j. cosy i Vcosy 
Vr 2 ' r 

(3) 

er = V siny 


(4) 


The perturbation parameter, e, which has been artificially inserted, is nominally one. It is assumed 
that the atmosphere is stationary and that the thrust vector lies in the vehicle's plane of symmetry. 
In (1), mass specific energy, E, is employed as a state variable in place of velocity, V, where 

E = V 2 /2 - p/r (5) 

The reference point for zero gravitational potential is taken at a radial distance approaching infinity. 
The symbol V is to be taken as 

V = [2(E + p/r)]l/2 (6) 

everywhere it appears in this analysis. Position and heading dynamics are decoupled from (1—4) 
by the assumption of a non-rotating Earth and are not of interest at present. 

Drag is assumed to have a conventional parabolic form 

D = qsC Do + KL 2 /qs where q = pV 2 /2 (7) 
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The assumed functional dependence for Cd 0 » the zero lift drag coefficient, and K, the coefficient 
of the induced drag component, are: 


C Do = C Do (M) K = K(M,a) 

Lift is given by 

L = qsCL(r,E,a,5 e ) = qsCLa(a - <*zl) 


( 8 ) 


(9) 


The lift curve slope, CL a . and the angle of attack for zero lift, azL > are assumed to be Mach 
number dependent. 


A multi-mode propulsion system composed of n different engine types (i.e. cycles) is assumed. 
Net thrust is given by 

T = [ F c 2 + F s 2 ] 1/2 (10) 


where Fc represents the component of net thrust along the velocity vector and F s represents the 
component of net thrust normal to the velocity vector, i.e. in the lift direction. These components 
are depicted in Figure 1 and given by: 

P A 

Fc = X Uj 1j cos (a + £Tj ) + X Ti cos (a + ) (11) 

j=l i=P+ 1 

P n 

Fs=XTljT j sin(a + er i )+ X Tisin(a + e Ti ) (12) 

j=l i=p+l 


Each of the n engine cycles (i.e. turbojet, ramjet, scramjet, rocket, etc.) is controlled by variation 
of the fuel flow rate in direct proportion to command. Of the total number of engine types to be 
considered, p are assumed to exhibit a linear relation between fuel flow rate and thrust generation. 
Each engine of this type shall be controlled by varying its throttle setting, T]j. This assumption is 
typically employed for rocket engines. For the remaining n-p engine cycles, the relation between 
fuel flow rate and thrust generation is assumed nonlinear. Control of each engine of this type shall 
be effected by variation in its fuel equivalence ratio, (pi* This behavior is typical of air-breathing 
cycles. The subscripted symbol T k (k = 1 to n) in (1 1) and (12) represents the net thrust generated 
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by any number of independent engines employing a particular cycle k. The symbol ejk denotes the 
angle between T k and the body longitudinal axis (see Figure 1). Note that in general. 


Tj = Tj(r,E,a) 

m 

H- » 

1 t 

j = 1 top 

(13) 

Tj = T; (r,E,cpi,a) 

CPi€ [0,1] 

i = p + 1 to n 

(14) 

Ej ~ £j (r,E,a) 


k = 1 to n 

(15) 


The total fuel flow rate, f, is given by 

P v 1 

f= / ,Tij Cj(r,E,a)Tj(r,E,a)+ ^ tPi fmaxi ( r ,E,(x) (16) 

j=l i=p+l 

where Cj represents the thrust specific fuel consumption for engine type j and fmax; represents the 
product of thrust specific fuel consumption and thrust at a stoichiometric fuel-to-air ratio for engine 
type i. For convenience all of the engine throttle controls are collected into a single vector as 
follows. 


rc T = [Til, Tl2, ... , Tip, <p p+ i, ... , <p n ] ( 17 > 

The control variables are angle of attack, a, bank angle, <7, the fuel equivalence ratios, <pi, for 
engine types 1 through n, and engine throttle settings, T|j, for engine types n+1 through p. The 
objective is to minimize the fuel consumed in gaining energy, with the performance index given 
by. 


J = -m(tf) (18) 

The final time, tf , is free. Minimization of (18) is to be carried out subject to maximum dynamic 
pressure and maximum aerodynamic heating rate inequality constraints and acceleration limits 
defined by 


Ci (r, E) — q qmax — 0 


(19) 
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C 2 (r, E, a) = Q-Qmax^O 


( 20 ) 


C 3 (r, E, m, cl, k) — nj ni max ^ 0 (21) 

C 4 (r, E, m, a, k ) = n 2 - n 2max < 0 (22) 

The symbols ni and n 2 represent the accelerations in g's along and normal to the velocity vector 
(i.e. in the lift direction), respectively. 


6 



SECTION 3 


Singular Perturbation Analysis 


3.1 Reduced Solution. 

Lift as a Control. We first consider a simplified problem in which flight is constrained to a 
vertical plane, the thrust vector is aligned with the velocity vector, and thrust production is 
assumed independent of vehicle angle of attack: 

a = 0, F S = 0, F C = T (23) 

Furthermore, we consider only that portion of the trajectory in the hypersonic regime. In this 
regime we need only consider a dual-mode propulsion system (i.e. n = 2). The system consists of 
a bank of scramjet engine modules assumed to operate continuously and a rocket engine that can be 
throttled as desired. The constraint (21), which can lead to the requirement for intermediate values 
of throttle setting, will be ignored. In this simplified setting the total fuel flow rate and net thrust 
can be represented as 

T = T s (r,E) + q T r (r) ; q e [0,1] (24) 


f — c s (r,E) T s + q c r (r) T r (25) 

where thrust specific fuel consumption is represented by c„ for the scramjet and c r for the rocket. 
Under these assumptions the governing equations of motion can be written as. 


E = 


V(T — D) 
m 


(26) 


m = -f (r,E,q) 


(27) 


* L U cosy V cos Y 

67 = mV 77T~ + r 
V r 


(28) 


er = V sinY 


(29) 
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The control variables are now rocket engine throttle, T|, and vehicle lift, L. The objective remains 
to minimize the fuel consumed in gaining energy. 

Setting 8 = 0 in (26-29) reduces the order of the dynamic system to two and results in what is 
conventionally referred to as the energy state approximation. That is, altitude and flight path angle 
dynamics are assumed fast in comparison to energy and mass dynamics, and altitude now takes on 
the role of a control variable 9 . The differential equations (28) and (29) are reduced to algebraic 
equations which yield the following relations: 

Yo = 0 (30) 


L 0 = m[(|i/r 2 ) - (V 2 /r)] 


(31) 


The subscript zero denotes reduced solution values and is omitted below where not deemed 
necessary for clarity. The reduced solution Hamiltonian is given by 


where 


H 0 = Xe E + A, m rti + constraints = 0 


1-0 


(32) 

(33) 


Satisfaction of the minimum principle with respect to altitude, h, is equivalent to the following 
operation (see Appendix B of reference 10), 


h 0 * = arg max h [V(T - D)/f] 


E = constant 
T > D q < qmax 

T|=T|* Q< Qmax 


(34) 


Consideration of the constraints (19, 20) simply limits the search space over which the 
maximization of (34) takes place. The superscript asterisk denotes an optimal value of control. 
This operation yields an optimal altitude program as a function of vehicle energy and mass. Note 
that T| appears linearly in the Hamiltonian resulting in a bang-bang control solution for rocket 
throttle setting. A switching condition, S, results from the evaluation of 3 Ho/<)t| and is given by, 

S=X E (V/m)-X m c r f (35) 
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Using (32) to eliminate X. E in (35) and taking into account the sign of X m + yields the following 
analytic switching condition: 


r| = 0 if [(Cj - C s )/Cr]T s > D 
T1 = 1 if [(Cr - CsVcJTs < D 


(36) 


Intermediate values of rocket throttle setting are not optimal. This fact is revealed by examination 
of the matrix H uu , which is required to be at least positive semidefinite along an optimizing singular 
arc. For convenience, V, rather than h is taken as the control-like variable so that u T = [ V,T) ]. The 
determinant of H uu , which is symmetric, must be greater than or equal to zero for positive 
semidefiniteness. However, it can be shown that 


det H m =-{H Vn | 2 


(37) 


which is negative for Hv n * 0, which is generally the case. 

It can happen that the velocity set is not convex in a region of interest, and, in the absence of 
convexity, one can not guarantee that an optimal control exists. Thus the possibility of a chattering 
control solution should be examined. Conclusions regarding this matter are model dependent and 
are discussed in reference 16. It is sufficient to say here that no chattering solutions for rocket 
throttle setting were found for the vehicle models examined. 

The reduced solution costates are determined as 12 : 


o 

II 

J 

^m 0 = - m(t f )/m 

(38) 

ho = ^E 0 [2KV 2 W(qs)] 

^-E 0 = km 0 (fm/[V(T-D)]) 

(39) 


t The "influence function", X m , represents the variation in the performance index, J, with 
respect to mass 11 . Since J = - m(tf), h cannot change sign (i.e. it is not possible for a reduction 
in vehicle mass as fuel is expended along the climb path to increase the final mass of the vehicle). 
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Angle of Attack as a Control . Consider now the full model complexity formulated in Section 
II with the exception that flight remains constrained to a vertical plane. That is, consider flight over 
the entire Mach range, including the subsonic and supersonic regimes. Assume a multi-mode 
propulsion system consisting of turbojet, ramjet, scramjet, and rocket cycles (i.e. n = 4). Allow for 
a component of net thrust normal to the velocity vector and consider the possibility that the 
performance of one or more of the air-breathing engine cycles is dependent on vehicle angle of 
attack. Consider also the constraint on axial acceleration given by (21). The method of solution 
proceeds as before. 

Setting 8 = 0 in (1-4) reduces the differential equations (3) and (4) to algebraic relations: 


o 

II 

(40) 

L 0 = m(p/r 2 -V 2 /r)-F s 

(41) 


The control oto is eliminated via (41). That is, given values of r and E, cq, is iteratively determined 
using (41) while enforcing trimjhrough elevon deflection, 8 e . More concisely, 

a 0 (r,E) = {a 0 : L(r,E,a,5 e ) - L 0 (r,E,7t,a) = 0} (42) 

The reduced solution Hamiltonian is again given by (32). But since drag, given by (7), is 
dependent on L 0 2 , which in turn depends on engine controls through Fs, as given in (12), the <pi 
and the T|j both enter nonlinearly in the Hamiltonian. Satisfaction of the minimum principle with 
respect to h and 7t is equivalent to the following operation, 


h 0 *, Jt* = arg max h % [V(Fc - D)/f] 


E “ constant 
q < qmax 

n l - n lmax 


FoD 
Q< Qmax 

n 2 ^ n 2max 


(43) 


This operation yields both an optimal altitude program as a function of vehicle energy and mass 
and the corresponding optimal engine controls. 


If we neglect the dependence of reduced solution drag on the sine component of net thrust, Fs, 
then the rij enter linearly in Hq. In such case we have bang-bang solutions for the T|j with possible 
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singular arcs along which intermediate throttle settings may be optimal. The switching functions 
are determined as before from dHo/cfrjj, 

Sj = [f(7t) cos(a + eTj)/(Fc(rc) - D)] - Cj j = 1 +n to p (44) 

Throttle settings are then governed by the following relations: 

T|j = T"lmin if Sj < 0 

T) j singular if Sj = 0 for finite time (45) 

l"lj = ^lmax if Sj > 0 

The Sj are dependent on the r^, k * j. Thus an iterative scheme is required to arrive at the optimal 
combination of throttle settings if j > 1. 

Thus far in the analysis it has been assumed that each engine cycle can be independently 
controlled. Since much of the captured mass flow and some or all of the engine hardware will be 
shared by the various engine cycles employed, it is perhaps more useful to consider operation of 
the various air-breathing cycles as mutually exclusive. In reality, dual combustion over a finite 
range of Mach number will be required to smoothly transition from subsonic to supersonic 
combustion 13 . One can view the case of mutually exclusive engine cycles as a problem in which 
the system equations are discontinuous at cycle transition points along the trajectory. Following the 
terminology of reference 1 1 , suppose that one set of system equations. 


X = f (1) (x,u,t) (46) 

applies for t < ti, where ti is free, and another set of system equations applies for t > q, namely, 

x = f< 2 >(x,u,t) (47) 

Here x and u denote general state and control vectors, respectively. It is necessary for optimality 
that 

H(D (ti-) = H(2) (t!+) (48) 

The condition (48) can be used to determine the optimal point of transition from one set of system 
equations to another. In this case f (1) and f (2) differ only by the thrust produced by the particular 
engine cycle being employed and by the associated difference in fuel consumption. Satisfaction of 
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(48) can be reduced to the following equality where the condition Hq = 0 has been employed to 
eliminate costate dependence 

(Ti cos(a + £tj) - D)/c,Ti = (Tj cos(a + ejj) - D)/cjTj (49) 


This result is, in fact, obvious from examination of (43). That is to say, points at which a change 
in engine cycle can occur require that the function to be maximized be equal for either choice of the 
propulsion cycle. When the functional evaluations are not equivalent, one or the other is greater 
and dictates the optimal choice of cycle. 


Bank Angle as a Control. It is reasonable to assume that the performance of the proposed 
scramjet engines will be sensitive to vehicle angle of attack. Furthermore, it is quite likely that 
thrust production will depend on angle of attack in a nonlinear way. Given that this is true, any 
particular engine installation will exhibit an angle of attack for which engine performance is best. 
This angle of attack for best engine performance, call it the design angle of attack, may in turn vary 
with Mach number 14 . If such nonlinear behavior is assumed and the optimal flight path is 
constructed using (43), one finds that, since fuel optimization is very sensitive to engine 
performance, the optimal trajectory tends to remain on a contour along which the design angle of 
attack is maintained. It can happen, however, that overall performance is improved if the design 
angle of attack is maintained while flying at lower altitudes, and hence at higher values of dynamic 
pressure. Of course, maintaining the design angle of attack at a higher dynamic pressure generates 
additional lift which causes the vehicle to immediately climb above the desired flight path. Thus, in 
order to fly along the optimal path, the extra lift associated with maintaining the design angle of 
attack must be "dumped." One procedure for accomplishing this task is to roll back and forth in 
such a way that, on average, the component of lift in the vertical direction is reduced to that 
required to maintain the optimal climb rate. It may in fact be more practical to appropriately offset 
the initial vehicle heading and to then execute a single coordinated turn that accomplishes the same 
objective. With bank angle thus introduced as an additional control, satisfaction of the minimum 
principle with respect to h, 7t, and a is equivalent to the following operation, 


ho i Tt * Oto 


= arg max h r a [V(Fc - D)/f] 


E = constant 
q < qmax 

n l - n lmax 


FoD 
Q< Qmax 

*2 - n 2max 


(50) 
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where bank angle, a, is determined so that 

L cos a = L 0 


(51) 


3.2 Boundary Layer Analysis. 

The unconstrained boundary layer solution associated with (1—4) is obtained by introducing 
the time transformation x = t/e and again setting e = 0. That is, energy and mass are held constant 
while altitude and flight path angle dynamics are examined on a stretched time scale. The resulting 
necessary conditions for optimality yield an optimal feedback guidance law for lift control which 
depends on the unknown costate Xy (see Section 5.2.2 of reference 10). In the absence of a state 
inequality constraint such as (19), a suitable approximation to X^can be obtained by linearizing the 
boundary layer necessary conditions about the reduced solution 10 ’ 15 . However, this procedure is 
not applicable when the reduced solution lies on a state constraint boundary. This problem is 
discussed in detail in Appendix A. The boundary layer control solutions for engine throttle settings 
are similar to those of the reduced solution. The T\j enter the Hamiltonian linearly, but the 
switching conditions that govern their behavior are also dependent on the unknown costate Xy. 
This dependence drops out of the switching conditions if the sine component of thrust, Fs, is 
neglected. 


Feedback Linearization - Lift as a Control. As an alternative approach to handling the 
control of altitude and flight path dynamics, a nonlinear transformation technique is employed as 
follows 10 . Consider the boundary layer altitude and flight path angle dynamics given in (32) and 
(33) on a transformed time scale x = t/e. Note that we have system equations in block triangular 
form. To proceed we take successive total time derivatives of r until explicit dependence on the 
control appears. The prime notation denotes differentiation with respect to x. 

r" = (Lcosy)/m + (V 2 cos 2 y)/r - (|J./r 2 ) (46) 

The control, L, appears in the second time derivative and we define U, the pseudo control, as 


U = r" 


(47) 


It is desired that U be determined as follows 
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U=K p (r 0 -r) + K d (r„-r') 


(48) 


where r 0 denotes the reduced solution radius at the current energy level and the time derivative of 
r 0 denotes the climb rate required to stay on the reduced solution as energy is gained. This climb 
rate can be estimated by defining an appropriate increment in energy, evaluating the reduced 
solution at this higher energy level, and then estimating the required climb rate using a forwards 
difference. 

The inverse transformation is defined by solving for L in (47) using (46) and (48), 

L = (U + (p/r 2 ) - [(V^cos 2 '/] } (m/cosy) (49) 

This lift control solution is constrained directly by (21). Note that as r and y approach their 
reduced solution values, (49) approaches the reduced solution value of lift given by (35). A block 
diagram depicting the conceptual implementation of the nonlinear transformation technique to yield 
the controller defined by (49) is presented in references 10, 15 and 16. The corresponding closed 
loop transfer function is 


G(s) = (K d s + K p )/(s 2 + KdS + K p ) (50) 

where the gains K p and Kd for the second order system can be written in terms of the damping 
ratio, £, and natural frequency, tOn, as 

K p = C0n 2 Kd = 2^COn (51) 

The performance of this controller can be dictated by selecting the values of K p and Kd to yield the 
desired dynamic response. This lift control solution applies equally well to the unconstrained or the 
inequality constrained case. 


Feedback Linearization - Angle of Attack as a Control. Direct extension of the lift 
control solution presented above to include the angle of attack effects included in (1-4) results in 
the following feedback control law. 
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The pseudo-control, U, is defined as before where again K p , proportional gain, and K<j, rate gain, 
are selected to yield the desired controller performance. Optimal lift, which is directly constrained 
by (21), is then determined by, 

L* = qsC^ (a* - cxzl) (64) 
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SECTION 4 
Vehicle Models 

Four different vehicle models were employed to generate the numerical results presented in the 
next section. The first, referred to as Model 1, is based on a hypersonic research vehicle concept 
studied by NASA in the 1970's and is powered by a combination of scramjet and rocket 
propulsion 10 . This model is useful only in the hypersonic regime. Models 2 and 3 are based on a 
"Generic Hypersonic Aerodynamic Model Example," or GHAME, developed more recently by 
NASA 17 . A nominal configuration of 233.4 feet total length and 300,000 pounds gross take-off 
weight was assumed. The trimmed aerodynamic characteristics were taken directly from the 
GHAME documentation. For Model 2 the largely empirical GHAME I aerodynamic data set was 
employed. For Model 3 the numerically generated GHAME II aerodynamic data set was employed. 
Both sets extend from take-off to orbital velocities. Thrust for both Models 2 and 3 is provided by 
a multi-mode propulsion system composed of turbojet, ramjet, scramjet, and rocket engines. The 
airbreathing propulsive characteristics for this model were adopted from reference 18. A rocket, 
sized for orbital insertion (roughly 15,000 lbs. of thrust in vacuum), is assumed available over the 
entire Mach range 10 . This system corresponds to the case p = n = 4 in (1 1), (12), and (16). As a 
result, the switching conditions given by (45) can be used to determine all of the cycle transition 
points. Figure 2 presents the adopted variation in fuel specific impulse with Mach number for the 
various engine cycles. The various engines where sized by trial and error and do not represent an 
optimal configuration. The generation of scramjet thrust due to mass ejection when operating above 
a stoichiometric fuel-to-air ratio is not modeled. Thrust induced pitching moments, which can be 
significant 16 , were not considered when trimming the aircraft. A fourth model was constructed by 
combining an aerodynamic data set provided by the NASA Langley Research Center (referred to as 
the "Langley Accelerator") with the propulsive data described above. Additional details regarding 
these models are available in references 10 and 19. 

A simple model for convective heating rate per unit area, Q, was adopted from reference 20, 

Q = (4.919 E-08) p 0 - 5 V 3 0 (52) 

Equation (52) gives Q in Watts/cm 2 given density in kg/m 3 and velocity in m/sec and corresponds 
to equilibrium conditions on the surface of a sphere or wing leading edge 10 cm in radius and 
cooled by reradiation alone. For reference, a contour of Q = 800 in the altitude-velocity plane 
corresponds roughly to a contour along which skin temperature remains at approximately 2000° F 
three feet aft of the leading edge assuming laminar flow 21 . 
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SECTION 5 

Numerical Results 


Reduced solution trajectories were generated by carrying out the maximization process 
indicated in (34), (43), and (50) over the energy range from take-off to orbit. Numerous results are 
available in reference 19. Only a few representative plots are presented here. Figure 3 depicts 
reduced solution trajectories for Model 1 in the altitude/velocity plane. Dynamic pressure is limited 
to 2000 psf while maximum allowable heating rate is varied. The trajectories follow the dynamic 
pressure constraint boundary until the specified contour of maximum heating rate is encountered. 
The path then follows the constant heating rate contour until reaching the trajectory for which no 
heating rate constraint was enforced. At this point the heating constraint becomes inactive and the 
trajectory rejoins the unconstrained climb path. The mechanism causing the altitude discontinuity at 
a velocity of 22,000 ft./sec. is similar to that which has been noted in the transonic region for 
supersonic fighter aircraft 15 . Included in Figure 3 is the rocket switching surface, i.e. the contour 
along which the switching function (39) remains zero. At altitudes below this contour the optimal 
rocket throttle setting is zero whereas above the contour the optimal throttle setting is one. The 
performance penalty paid in enforcing the heating constraint is presented in Figure 3 as time 
required and percent gross weight consumed to achieve an orbital energy level. This performance 
penalty must be weighed against the complications of using active cooling, the weight of heat 
shielding, and various other factors in the vehicle design process. 

Figure 4 presents the reduced solution climb path in the altitude-velocity plane for Vehicle 
Model 4. The dynamic pressure constraint is again enforced, an aerodynamic heating constraint is 
not, and a limit on axial acceleration is introduced. The trajectories vary predictably as the 
magnitude of the acceleration constraint is changed. Some throttling of the engines is employed but 
in general the vehicle prefers instead to climb to reduce the excess thrust available. Note that the 
altitude discontinuity present in Figure 3 does not occur for this vehicle, but is in fact implicit in 
matching a terminal altitude condition at orbital velocity. The horizontal bars at the top of the figure 
indicate the velocity range over which the operation of each engine cycle was deemed optimal, 
including regions of cycle overlap. The rocket operation was not optimal during atmospheric flight 
for this case. 

Computational investigations of the sensitivity of scramjet performance to changes in angle of 
attack predict highly nonlinear behavior 14 . Figure 5 presents a scramjet thrust scaling factor 
employed to model this effect. This figure is based on a liberal interpretation of the computational 
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results presented in reference 14. This curve is shifted with respect to the horizontal axis in order 
to represent inlet designs which favor maximum engine performance at an angle of attack other 
than zero. 

Figure 6 presents variations in the dynamic pressure constrained reduced solution trajectories 
for Model 3 when thrust variation with angle of attack for both the ramjets and scramjets, as 
depicted in Figure 5, is included. The angle of attack for best engine performance is varied over the 
range from 0.5 to 3.0 degrees. When only scramjets are operating, the vehicle tends to prefer a 
path along which the "design" angle of attack is strictly maintained. The performance penalty paid 
for a change in the design angle of attack is modest however, since in this Mach region, the 
acceleration capability of the vehicle is high. When the thrust scaling factor of Figure 5 is assumed 
Mach dependent in accordance with the results of reference 14, a much greater variation in the 
trajectory is experienced 19 . 

The peak in the trajectories at approximately 3000 ft./sec. in Figure 6 is due to turbojet shut 
down. This peak is significantly reduced when the turbojet inlet area is increased, indicating that 
the climb away from the dynamic pressure constraint boundary is due to the decreasing level of 
thrust available from the turbojet as the Mach number increases. With an increase in altitude comes 
a reduction in vehicle drag, but the turbojet switching surface is encountered at an altitude of 
approximately 75,000 ft. and the turbojet shut down. The SCRAMJET almost immediately 
switches on, and with a much greater magnitude of thrust, can sufficiently overcome vehicle drag, 
even at a higher dynamic pressure. Thus the trajectory returns to the dynamic pressure constraint 
boundary. Note that the ramjet is turned on at a very low Mach number (i.e. M = 0.81) even 
though it is extremely inefficient in this speed range (see Figure 2). This behavior has been noted 
by past researchers and is due to the presence of a "pinch point" (i.e. a point of minimum thrust 
minus drag) in the transonic region. The size of the ramjet was selected without regard to its 
weight. However, optimization of the vehicle configuration must take into account the mass of 
each engine and the mass of the required engine cowling. Results obtained indicate that the optimal 
trajectory for such an optimized configuration may prefer the use of rocket (rather than ramjet) 
thrust to augment turbojet thrust at the transonic pinch point. 

As stated above, cycle operations are represented in Figure 6 by horizontal bars. The transition 
points were very nearly the same for thrust independent of or dependent on angle of attack. The 
overlap in air-breathing cycles is desirable to provide smooth cycle transitions. For Model 3, 
turbojet sizing requires about 25 sq. ft. of inlet area, whereas the total number of ramjet modules 
selected require 200 sq. ft. of inlet area. Thus it should be possible to start the majority of the 
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ramjet engines in a sequence that avoids excessive accelerations while maintaining turbojet thrust. 
Once the velocity for turbojet shut-down is reached, turbojet airflow can be diverted to the 
remaining ramjet modules. The number of required scramjet modules is likewise larger than the 
number of required ramjet modules, and a similar argument for cycle overlap can be made. Of 
course, the actual system will no doubt share much of the engine hardware amongst the various 
cycles employed in addition to sharing the captured mass flow. Thus the actual optimization of 
engine transitions will be more complex. 

Modulation of the vertical component of lift via bank angle variation was also evaluated for 
Model 3. Carrying out the maximization process indicated in (50) alters the trajectories presented in 
Figure 6 only slightly. The changes correspond to those portions of the trajectory where cto < 
“design. As such, only a very modest gain in performance was achieved. However, if design 
constraints force the scramjet design angle of attack to differ significantly from the angle of attack 
for zero lift, much greater savings can be obtained. 

Figure 7 depicts the reduced solution trajectory for Vehicle Model 2 in the altitude-velocity 
plane. A maximum allowable dynamic pressure of 2000 psf is the only constraint enforced. The 
dashed line labeled 1 represents the fuel-optimal climb path when scramjet performance is assumed 
independent of vehicle angle of attack. The percent of take-off gross weight consumed in attaining 
orbital energy is 61. The solid line label 2 represents the fuel-optimal climb path when scramjet 
performance is assumed to vary with angle of attack according to Figure 5, with optimum engine 
performance assumed to occur at an angle of attack of 3 degrees across the Mach range. In this 
case the trajectory tends to remain on the dynamic pressure constraint boundaiy for the majority of 
the flight The percent of take-off gross weight consumed in attaining orbital energy in this case is 
68.2. The weight penalty of 7.2 percent of the take-off gross weight most likely exceeds the 
payload capability of the vehicle. This comparison indicates the critical need to accurately model 
the many interactions present among disciplines. 

Figure 8 depicts the altitude time history for simulated flight of Model 3 using the lift control 
law derived via feedback linearization to track the corresponding reduced solution. The ramjet 
cycle was eliminated and the trajectory is subjected to the following constraints: dynamic pressure 
<» 2000 psf, reference heating rate < 400 Watts/cm 2 , and axial acceleration < 3.0 g's. In general 
this vehicle preferred to climb in order to satisfy an axial acceleration limit rather than to throttle 
back the engines. The rapid climb at roughly 400 seconds is due to scramjet tum-on and this 
preferred behavior. The large overshoot just before 500 seconds is due to the inability of the 
vehicle to pull down as the altitude for which ni < 3 at full throttle is approached. This overshoot 
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can be reduced somewhat by careful gain scheduling and having the controller look ahead in 
energy. However, the requirement to fine tune the controller for each trajectory generated is not 
desirable for the intended applications of this algorithm. The problem has less to do with the 
controller than with the generation of the trajectory itself. 

Over the vast majority of the trajectory the flight path angle is small and the flight path angle rate 
very modest so that (41) provides a good approximation to the actual lift required to follow the 
flight trajectory. However, when the scramjet is initially turned on at a relatively high value of 
dynamic pressure, the energy rate of the vehicle is greatly increased. The necessity of 
simultaneously climbing to avoid violating the dynamic pressure constraint boundary results in a 
large flight path angle rate. The time scale separation assumed in (1-4) is simply not appropriate 
over this small portion of the trajectory. A simple way to overcome this difficulty consists of 
estimating the flight path angle and time interval between energy levels, combining them to form an 
estimate of the flight path angle rate, and then inverting relation (3) to obtain the required lift. By 
restricting the accelerations normal to the flight path when constructing the reduced solution in this 
region, a feasible trajectory can always be obtained. 

Figure 9 depicts the reduced solution climb path for Model 3 again with a maximum dynamic 
pressure of 2000 psf, a maximum aerodynamic heating rate of 400 Watts/cm 2 , but with a 
maximum axial acceleration of lg to amplify the problem (y = 0 in the lift calculations). Also 
depicted is the modified trajectory when the method described above is implemented (y and dy/dt * 
0). The results in the altitude/velocity plane are quite dramatic over the speed range from 3,000 to 

12.000 ft./sec. The near vertical altitude transition at a velocity of approximately 3,000 ft./sec. is 
eliminated, as is the dive that followed. The arc which follows in the velocity range from 5,000 to 

13.000 ft./sec. corresponds to the region over which the axial acceleration limit is active. Less 
altitude change is commanded in this region; more throttle is used to reduce the axial acceleration 
instead. The remainder of the trajectory, the same for either case, constitutes flight along the 
heating constraint boundary. Despite the significant change in trajectory, only 200 additional 
pounds of fuel are consumed and the difference in time of flight is only about 60 seconds. These 
small differences are due to the fact that the velocity interval from 3,000 to 12,000 ft./sec. is 
traversed very rapidly in time, corresponding to only a small fraction of the total time of flight. 
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Figure 1 Partial Force and Moment Diagram 


Turbojet 



(jjs) asjndiu j aifpads pnj 


22 


Figure 2 Fuel specific impulse versus Mach number for the various propulsion cycles available to Vehicle Models 
2 through 4. 
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Figure 3 Reduced Solution for Vehicle Model Number 1 as the maximum aerodynamic heating rate constraint 
(shown in Watts/cm 2 ) is varied. 
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Figure 4 Reduced Solution for Vehicle Model Number 4; variation in the magnitude of the axial acceleration 
constraint is shown. 
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Figure 5 Assumed variation in scramjet thrust with angle of attack. 
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as the design angle of attack is varied. 
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Figure 7 Reduced Solution for Vehicle Model Number 2; the effect of thrust variation with angle of attack is shown. 
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Figure 9 Reduced Solution for Model 4; comparison of the trajectory when flight path angle and flight path angle 
rate are included in the calculation of lift with the trajectory when they are not included. 



SECTION 6 

Conclusions and Recommendations 


6.1 Conclusions 

This research effort has demonstrated the utility of singular perturbation methods in the study of 
single-stage-to-orbit airbreathing vehicles and, in particular, in the derivation of efficient algorithms 
for ascent trajectory optimization and optimization of engine cycle transitions. The analysis extends 
over the entire Mach range from take-off to orbit and accomodates a realistic nonlinear vehicle 
model and all pertinent trajectory constraints. A number of important modeling and analysis issues 
not treated in the early stages of this effort were identified and addressed during this reporting 
period. Reasonable assumptions regarding propulsion system characteristics were introduced that 
allow the optimal engine cycle transition points to be determined as a function of state using a 
simple iterative test. These switching conditions lead to significant computational savings during 
the optimization process. Functional dependence of scramjet thrust on vehicle angle of attack was 
shown to have a major impact on the nature of fuel-optimal ascent trajectories. Also, depending on 
the actual vehicle configuration and the characteristics of the engine inlets, roll maneuvers used to 
modulate the vertical component of lift were shown to sometimes improve the index of 
performance during ascent. Over those limited regions of flight where the energy state 
approximation was found to be poor, simple lift corrections that account for non zero flight path 
angle and flight path angle rate were introduced that significantly improve the trajectory generation 
methodology. 

6.2 Recommendations 

Future efforts should be directed towards enhancing the performance and applicability of the 
derived algorithm. Such efforts should include the development of detailed multi-disciplinary 
vehicle models, the incorporation of additional controls such as thrust vectoring, reaction jets, and 
variable geometry, optimal control of the total heat load on the vehicle, the study of three- 
dimensional maneuvers, including abort, an examination of robustness issues, and improvements 
in speed of operation. 

The demonstrated capability for rapid near-optimal trajectory generation has yet to be exploited 
in the development of efficient tools for fully integrated hypersonic vehicle design. Efforts to move 
in this direction should include tieing a parameter optimization algorithm around the trajectory 
optimization code that has been developed and the incorporation of this algorithm into a fully 
integrated control system design methodology. 
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Work of more general theoretical interest was also initiated during this reporting period. It was 
found that state-variable inequality constrained boundary layer systems are not well understood. 
Many of the characteristic features of such systems were identified. For instance, it was found that, 
when the reduced solution lies on a state-variable inequality constraint boundary, the boundary 
layer trajectories are of finite time in the stretched time scale. The possibility of costate 
discontinuites at the juncture between constrained and unconstrained arcs makes direct application 
of existing theory difficult. A transformation technique was identified that eliminates some of these 
difficulties, but at the cost of possibly increased system order and the introduction of singular arcs. 
Further research in this area is recommended. 

Continued work with the integrated aerodynamic/propulsion performance prediction program 
has resulted in a highly accurate and useful means both for providing the needed vehicle parameters 
in the present program and for more general transatmospheric flight performance calculations. The 
program is evolving into a completely interactive performance estimation package, which will make 
it possible to view effects of small configuration changes on any performance parameters. The 
user can view in animated graphical form the effect of desired vehicle configuration changes. 
These modifications can be entered graphically by moving defining points on the vehicle outlines 
or by means of shifting simulated “levers” built in to the computer program. For example wing 
incidence angle, twist, wing area, fin cant angle, etc. can be changed continuously with 
simultaneous graphical output showing the effect on selected performance parameters. We 
anticipate many applications for this analytical capability and will continue to improve upon it. 

6.3 Publications 

Four conference papers have now been published which discuss most of the results of this 
research effort 12 ’ 15 * 16 - 49 . A Ph.D. Dissertation that details the entire effort to date was published in 
Decemeber of 1989 19 . A full-length paper entitled "Rapid Near-Optimal Trajectory Generation for 
Single-Stage-to-Orbit Airbreathing Vehicles" has been submitted for publication in the AIAA 
Journal of Guidance, Control and Dynamics and a new paper is now being prepared for the 1990 
AIAA GN&C Conference on the issue of state contraints in singularly perturbed systems. 
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Appendix A 

State Inequality Constrained Boundary Layers 


Abstract 

The established necessary conditions for optimality in nonlinear control problems that involve 
state-variable inequality constraints are applied to a class of singularly perturbed systems. The 
distinguishing feature of this class of systems is a transformation of the state-variable inequality 
constraint, present in the full order problem, to a constraint involving states and controls in the 
reduced problem. It is of particular interest to construct the zeroth order boundary layer solution 
when the reduced solution lies on the constraint boundary. It is shown that, in general, the 
boundary layer problem is of finite time in the stretched time variable. A special case is identified 
in which the boundary layer time scale transformation results in an increase in state inequality 
constraint order. In this case, required smoothness properties possessed by the full order system 
may be lost, and the application of existing necessary conditions for singularly perturbed systems 
then becomes invalid. A Valentine transformation can be used to regain required smoothness, but 
at the price of introducing singular arcs and an increase in system order. Finally, the various 
system properties and characteristics described in the body of the appendix are illustrated with 
several simple examples. 


I. Introduction 

State inequality constraints are commonly encountered in the study of dynamical systems. The 
study of rigid body aircraft dynamics and control is certainly no exception. For instance, a 
maximum allowable value of dynamic pressure is usually prescribed for aircraft with supersonic 
capability. This limit is required to ensure that the vehicle's structural integrity is maintained and 
constitutes an inequality constraint on vehicle state. State inequality constraints have been studied 
extensively by researchers in the field of optimal control, and necessary conditions for optimality 
when functions of state are constrained have been obtained 1 * 3 . However, the construction of 
solutions via this set of conditions proves very difficult, and most practitioners rely on direct 
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approaches to optimization that employ penalty functions for satisfaction of state inequality 
constraints 4 . 

As discussed in the literature, the use of singular perturbation techniques in the study of aircraft 
trajectory optimization can, through order reduction, lead to both open and closed loop solutions 
that are computationally efficient. These methods can also be used to circumvent difficulties 
associated with enforcing a state inequality constraint in the reduced solution 5 . As an example 
consider the minimum time intercept problem of 6 . A near optimal feedback solution is obtained via 
singular perturbation theory that includes consideration of an inequality constraint on dynamic 
pressure. In the zeroth-order reduced solution, algebraic constraints are obtained when the 
perturbation parameter, e, which premultiplies the so called "fast" dynamic equations, is set to 
zero. These constraints can be used to eliminate the fast states (in this case altitude and flight path 
angle) from the reduced problem. One can choose, however, to retain one or more of the fast states 
and to eliminate instead some of the original control variables. The retained fast state variables are 
treated as new controls, and the original state constraint becomes a constraint involving both state 
and control in the reduced problem. In subsequent analysis of boundary layers, altitude resumes 
its status as a state variable, and dynamic pressure once again becomes a function of state alone. 
However, because the reduced solution for the example F-8 aircraft does not lie on the dynamic 
constraint boundary during ascent, the inequality constraint on dynamic pressure was not 
considered. Of note is the fact that modem supersonic fighter aircraft (such as the F-15) do ride the 
dynamic pressure constraint boundary during the ascent leg of the minimum time to intercept path. 

In addition to the example cited above, dynamic pressure bounds are encountered during fuel- 
optimal climb for supersonic transports 7 for rocket powered launch vehicles such as the U.S. 
space shuttle 8 , and for single-stage-to-orbit air-breathing launch vehicles 9 . If, as in applying 
singular perturbation methods in seeking a solution to any of these problems, the reduced solution 
climb path lies directly on the dynamic pressure constraint boundary for a portion of the flight, 
then it is necessary to consider boundary layer transitions onto the constrained arc. This problem, 
which proves quite perplexing, has received almost no attention in the literature. 

This appendix documents an initial investigation of the features of boundary layer transitions to 
state constrained arcs. Section II provides a brief review of first order necessary conditions derived 
for state-variable inequality constrained problems in optimal control. Section III discusses the 
optimal control of singularly perturbed systems subject to state-variable inequality constraints in 
general, and in particular examines the features of state inequality constrained boundary layers 
when the reduced solution lies on the constraint boundary. Section IV provides several simple 
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examples which illustrate the problem features discussed in earlier sections. Section V completes 
the appendix by providing some concluding remarks. 

II. Constrained Problems in Optimal Control 

The introduction of a state inequality constraint of the form 

S(x,t) < 0 (A-l) 

can lead to considerable difficulty when attempting to obtain an optimal control solution. One 
approach to incorporating state inequality constraints into necessary conditions for optimality 
consists of constructing successive total time derivatives of S until explicit dependence on the 
control appears 11 . If p time derivatives are required then (1) is referred to as a p th order state 
variable inequality constraint. The function SP(z,u,t)=0 is then adjoined to the Hamiltonian as a 
constraint to be enforced when S=0. This approach introduces the following additional tangency 
conditions at the point of entry to a constrained arc 


N(z,t) 


' S(z,t) 

sVz,t) 



Ls p_ 1 (z,t)J 


(A.2) 


These same tangency conditions also apply at a point where the path leaves the constraint 
boundary. The equations (2) constitute a set of interior boundary conditions that must be met at 
each juncture between a constrained and unconstrained arc. Unfortunately, in order to satisfy these 
interior boundary conditions one must allow for the possibility of discontinuities in the costate 
variables at the junctures. An alternative set of necessary conditions can be obtained by adjoining 
the constraint function, rather than its p 01 derivative, to the Hamiltonian and then employing a 
separating hyperplane theorem^. These conditions prove simpler and sharper than those of 
reference 1 1 however the possibility of discontinuous costates is still present. The gap between 
the necessary conditions of references 11 and 22 is defined in reference 23. A third alternative 
involves enploying a transformation technique in which a slack variable is used to transform the 
state inequality constrained problem into an unconstrained problem of higher dimension 30,31 . The 
work associated with the derivation of these first order necessary conditions is detailed in 
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references 32-38. Second order necessary conditions for optimality in the presence of state 
constraints have also been derived 3 ^ as have conditions for various special cases 40 . Work 
continues in the area of state constrained optimization as evidenced by the approach of Lowen 41 . 


III. Optimization of Singularly Perturbed Systems 
Subject to State Inequality Constraints 

Consider the system of singularly perturbed nonlinear differential equations: 


dx/dt = f(x,y,u,t) 

(A. 3) 

e dy/dt = g(x,y,u,t) 

(A. 4) 

with an index of performance of the form 


t, 

J = (p[z(t f ),t f ] + J[ L[z(t), u(t), t]dt 

(A.5) 


where x and f are of dimension n, y and g are of dimension m, x(t 0 ) and y(to) are given, e is a 
small parameter, to < t < tf, and the control u(t) is of dimension p. Zero order necessary conditions 
for optimality of the associated reduced and boundary layer problems in the absence of state 
constraints are readily available 42 . However, the following restrictions apply: f, g, df/dx, df/dy, 
dg/dx, and dg/dy must be continuous and u must be piecewise continuous. Because of these 
restrictions on smoothness a direct extension of the necessary conditions of for state constrained 
problems to include singularly perturbed systems is not possible. This is due to the previously 
mentioned fact that discontinuities in the costates can occur at the junctures between constrained 
and unconstrained arcs. Alternately, the state inequality constrained singularly perturbed problem 
of interest can be converted into an unconstrained singularly perturbed problem of higher 
dimension by introduction of a slack variable 30 . This approach does eliminate the problem of 
discontinuous costates. However, the state constrained arc is replaced by a singular arc and the 
prospect of increased system dimension is unwelcome given the basic tenet of seeking order 
reduction. 

Consider the flight dynamics problem detailed in the main body of this report. An inequality 
constraint on dynamic pressure of the form 
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S(h,V) = q-q max <0 


(A. 6) 


is to be enforced where 

q = pV 2 /2 (A.7) 

The symbol p represents atmospheric density and V represents the flight velocity. Given the 
equations of motion (1-4) expressed in the main body of this report, the first time derivative of 
(A.6) can be expressed as, 

dS/dt = dq/dt = [V 3 0p/3r)/2 - ppV/r 2 ] siny + pV(T - D)/m (A.8) 

Recall that the symbol T represents thrust, D, aerodynamic drag, m, vehicle mass, r, radial 
distance from the center of the Earth, y, flight path angle, and 1 1, the gravitational constant for the 

Earth. Assume, as is typically done, that atmospheric drag can be represented as follows, 

D = qsCcx, + KL 2 /qs (A.9) 

where s represents an aerodynamic reference area, C Do ,the zero lift drag coefficient, and K, the 
coefficient of the induced drag component. Note that the drag is explicitly dependent on the lift, L, 
which is treated as a control. In addition, the relation for thrust, T, is usually explicitly dependent 
on the engine throttle control. These controls appear explicitly in the first time derivative of the 
constraint function, (A.8), and it is thus classified as a first order state inequality constraint (i.e. p 
= 1). It is shown in reference 22 that when the constraint function is adjoined directly to the 
Hamiltonian and p = 1, no jumps in the costates will occur at the entry of an unconstrained arc onto 
a constrained arc. In this case the smoothness properties required by general singular perturbation 
theory are not violated and we may proceed with the application of singular perturbation methods 
with confidence. 

The state inequality constraint on dynamic pressure is conveniently reduced to a state and 
control constraint function in construction of the reduced solution (i.e p becomes zero). This 
occurs because altitude, a state variable in the full order problem, becomes a control variable in the 
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reduced solution. Enforcing this state/control constraint in the reduced solution is a trivial matter. 
Now consider the construction of the zeroth order initial boundary layer. Assuming that energy is 
characterized as a slow state, application of the time stretching transformation x = t/e and again 
setting e to zero yields an energy rate of zero. As such, the last term in (A. 8) is no longer present. 
Instead, the first time derivative of the constraint function is given by. 


dS/dx = dq/dx = [Y 3 (ap/3r)/2 - ppV/r 2 ] siny (A. 10) 


With the last term of (A.8) absent, control dependence does not explicitly appear in (A. 10). That 
is, the classification of the constraint function is altered following the time scale transformation 
when £ is set to zero. Taking the time x derivative of (A. 10), (i.e. forming the 2 nc * time x derivative 
of S), yields a term containing the time derivative of the flight path angle, y. The expression for 
flight path angle rate is as follows. 


L p cosy t V cosy 

Y “mV ~ 2 + r 
Vr 


(A.l 1) 


which is explicitly dependent on the lift control. Thus, in the boundary layer, the inequality 
constraint on dynamic pressure is 2 n( ^ order. Unfortunately, there is no guarantee that the costates 
are continuous for this case as there was for the case p = 1. Note that this type of behavior, in 
which the inequality constraint function order varies, is not present in all singularly perturbed 
problems with state inequality constraints, just for a certain class of them. For instance, if the 
constraint function is dependent on fast states alone, this variation does not occur. 


If jumps in the costates are in fact present at the juncture between an unconstrained and a 
constrained arc, the smoothness properties required are violated and the available necessary 
conditions for optimal control of a singularly perturbed system cannot be applied directly. In some 
cases it is not possible for the boundary layer dynamic system to "ride" the constraint boundary 
before reaching the reduced solution. In such case the boundary layer costates can, at most, be 
discontinuous at the juncture between the initial boundary layer and the reduced solution and only 
if the reduced solution at that point is on the constraint boundary. It is also possible that no such 
jumps will occur. 

If we proceed assuming that such jumps do not occur, then in most cases of interest the 
functional form of the boundary layer control solution in the presence of a state- variable inequality 
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constraint can be obtained. If this functional form involves an unknown costate, then use of the 
derived control solution is prevented unless the associated two point boundary value problem 
(TPBVP) is solved. To avoid solution of the TPBVP, as is desired in seeking a solution suitable 
for on-board, real-time implementation, an approximation for the unknown costate can be formed. 
If the costate history is continuous, the linearization technique of reference 43 can be applied to 
form the estimate. Unfortunately, for the case described above in which the reduced solution lies 
on the constraint boundary and the inequality constraint order is elevated to 2 by the boundary layer 
time scale transformation, purely imaginary roots result when the linearization technique is applied. 
Thus it is not possible to find a stablizing costate approximation given arbitrary initial states. In an 
unconstrained problem, the lack of an appropriate eigen-structure indicates that the problem does 
not exhibit the time scale properties assumed. Boundary layer transients do not exist for such 
cases; in fact there are fast oscillations that do not die out. The addition of an artificial cost term in 
formulation of the problem is suggested as an ad-hoc way to circumvent this difficulty. By proper 
choice of the weighting on can guarantee the proper structure of the linearized boundary layer 
system 44 . 

An interesting feature of the constrained boundary layer system described above is the presence 
of a finite costate rate at the juncture with the reduced solution. This behavior is illustrated in the 
sketch presented as Figure 1. In the figure, ti denotes the time at which an unconstrained arc joins 
a constrained arc (i.e. a juncture point). 
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Figure 1 Illustration of finite costate rate at juncture between boundary layer and reduced 
solution trajectories. 


Identification of the behavior illustrated in Figure 1 is based on the following construction. 


dH B L 

dy 


5H B l 


3y 

U = U* 

S<0 



u = constant 
S<0 



(A. 12) 


Note that, since the last term in (A. 12) is zero, this relation is equal to the negative of the costate 
rate. The following relation can also be constructed. 


43 


(A. 13) 


dH B L 

8y 


5H bl 


8y 

U = U* 

S< 0 


+ X < 


as 

3y 


U = u* 


unconstrained 


Now as x tends to tf~ we find that, 



u = constant 
unconstrained 




0 

^ * 


= -Xs 0 3S/3y | 

o 


at ti + this becomes Xs o 


(A. 14) 


Thus the fast costate rate, Xy', has a finite value at ti - and then jumps to zero at ti + . The boundary 
layer system no longer approaches the reduced solution asymptotically as in the unconstrained 
case. Instead, a finite time boundary layer is implied. Similar finite time boundary layer 
phenomenon have been discovered by the investigators of interior boundary layers and boundary 
layers that approach singular arcs 45 ’ 46 . Using this terminal value for the costate rate it is possible 
to show analytically for problems of interest that the linearized boundary layer system will always 
have purely imaginary roots when the reduced solution lies on the state-variable inequality 
constraint boundary. 

Consider the possibility of integrating the boundary system backwards in time from the reduced 
solution using the finite terminal value of costate rate to get started. Note that only a single extremal 
will be generated unless an additional free parameter is introduced into the problem. Since it should 
be possible to reach any set of initial conditions that do not violate the constraint, such a parameter 
surely exists. The only available parameter appears to be the magnitude of the possible costate 
jump at the juncture. This would imply that the costate history will be discontinuous at the juncture 
for all initial conditions that do not lie on the single extremal generated when it is assumed no jump 
occurs. 
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Assuming that first order necessary conditions for optimality can be obtained for which the 
requirement for smoothness can be relaxed (see for instance reference 47), one would again 
require a scheme for obtaining a stabilizing estimate of the unknown costates that appear in the 
optimal control law. However, to the author's knowledge, no directly applicable stability theory 
for finite time phenomenon is available for completing this task 48 . Alternately, one can consider 
transforming the constrained problem into an unconstrained problem, generally of higher order. 
The optimal trajectory of the transformed problem exhibits singular arcs which correspond, in the 
original constrained problem, to arcs which lie on the constraint boundary 30 . Because of this, the 
technique for costate approximation using the linearized boundary layer system, will fail. The 
reader is referred to the literature for a description of the work that has been done with regard to 
understanding control of singularly perturbed systems that include singular arcs 46 . 


IV. Examples 

Several simple examples are now presented which illustrates the application of the Valentine 
transformation technique described in reference 46 without the penalty of increased system 
order 31 . The phenomena of a finite time boundary layer is illustrated in example 2. 

Example 1 

Consider the following singularly perturbed dynamical system with initial conditions at zero. 


2 

Xj = x 2 -U 

xi(0) = 0 

(A. 15) 

ex 2 = x 3 

x 2 (0) = 0 

(A. 16) 

ex 3 = u 

x 3 (0) = 0 

(A. 17) 


The following 2 nd order (i.e. p = 2) state inequality constraint is to be enforced, 

S = x 2 - 1 < 0 (A. 18) 

The final value of x 1 is specified, the final time is free, and the performance index is given by, 
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/ ■ 

dt 

0 


(A. 19) 


Using Valentine's device, the inequality constraint (A. 18) is converted into an equality by the 
introduction of a "slack variable," a 30 . 

S + a 2 /2 = 0 (A. 20) 

Differentiating (A. 18) p times (p = 2) with respect to time, the following set of equations is 
obtained, 

X 3 /£ + acq/e = 0 where e da/dt = ai (A.2 1 ) 

u/e 2 + ai 2 /£ 2 + act 2/£ 2 = 0 where £ dai/dt = ct 2 (A.22) 


Using the transformations X 2 = 1 — <x 2 /2, X 3 = aai, and u = — <Xi 2 — aa 2 , (A. 15-17) become, 


a 2 / 2 f 


xi = 1 (ctj + aaJ 

(A.23) 

m 

EOL = a x 

(A.24) 

ea x * a 2 

(A.25) 

Reduced solution. By setting £ to zero, (A.24) and (A. 25) are reduced to the following, 

0 


otj = 0 

(A.26) 

a 2 = ° 

(A.27) 


The reduced solution Hamiltonian is given by, 
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>2 




1 _ T 

- ( a l 

+ aa 2 | 

+ + X a d 2 +1=0 

(A.28) 

Evaluation of first order necessary conditions for optimality results in the following. 





o 

0 


H = 0 & H a = 0 

=> 

K = ~ 

-1 , a = 0 

(A. 29) 

H Bi = 0 

=> 

0 

1 

(A. 30) 

Har 0 

=> 

0 

-1 

(A. 31) 


Boundary Layer Problem. Introducing the time scale transformation t = t/e and again setting £ 
to zero, the boundary layer dynamics are given by 

a = dj (A. 32) 

aj = a 2 (A.33) 

where the prime notation denotes differentiation with respect to the stretched time x. The boundary 
layer Hamiltonian is given by, 


h bl = 


a 


= + K + aa 2| + ^a a l + *-a«2= 0 


The costate dynamics are given by 



(A.34) 


(A.35) 
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(A.36) 





2 

ctj + aa 


a, 


Jt 


a 


And when the constraint is inactive, the control, 0 . 2 , is determined by the necessary condition, 


3H bl 

da 2 



2 1 

a x + aa 2 i 


2 

I a + 



(A. 37) 


Note that when the constraint is active, a is zero and the condition (A. 37) yields no direct control 
solution. Riding the constraint boundary corresponds to a singular arc in the transformed problem. 


Example 2 

Consider a simplification of Example 1, namely the singularly perturbed dynamical system, 

x j = x 2 ■ u^ (A. 38) 

ex 2 =u (A. 39) 

The inequality constraint (A. 18) becomes a 1 st order (i.e. p = 1) state inequality constraint, 

S = x 2 - 1 < 0 (A. 40) 

The final value of x is again specified, the final time is free, and the performance index is again 
given by. 


/ ► *i 

dt (A.41) 

0 

Using Valentine's device, the inequality constraint (A.40) is converted into an equality by the 
introduction of a "slack variable," a 30 . 

S + a 2 /2 = 0 (A. 42) 


48 



Differentiating (A.42) p times (p = 1) with respect to time, the following equation is obtained, 
u/e + aai/e = 0 where e da/dt = o.\ (A.43) 


Using the transformations X 2 = 1 — a 2 /2.and u = — acti, (A. 38-39) become. 



(A.44) 


ea = aj 


(A. 45) 


Reduced solution. By setting e to zero, (A.45) is reduced to the following. 


a 1 = 0 (A.46) 


The reduced solution Hamiltonian is given by, 


H = X, 


i - V - K) 


+ X a <x l +1=0 


(A.47) 


Evaluation of first order necessary conditions for optimality results in the following, 


H = 0 & H a = 0 => 

H a = 0 => 

a 2 

from which it is evident that 

x°(t) = t 


o o 

= -1 , a =0 

(A.48) 

0 

X = 0 

(A. 49) 

- 


(A.50) 
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(A. 51) 


*2(t) = 1 

Boundary Layer Problem. Introducing the time scale transformation x = t/e and again setting e 
to zero, the boundary layer dynamics are given by 


a=a i (A. 52) 

where the prime notation denotes differentiation with respect to the stretched time t. The boundary 
layer Hamiltonian is given by. 


h bl = 

The costate dynamics are given by 



(A. 53) 


A. a = -a(l + 2<xJ (A. 54) 

The condition that the partial derivative of the Hamiltonian with respect to the control be zero yields 
the following result. 


H a ,= 0 


a, = 



(A. 55) 


Substituting this result back into the condition that the Hamiltonian be zero, the following result is 
obtained. 


2 

H = 0 => 'k a = ± fl a. 


(A.56) 


Substituting (64) into (63) we find that the optimal value of the control, cti, is constant; namely. 


* 




(A. 57) 
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With this result it is possible to integrate equation (A.52) to obtain, 

a(t) = fl ± y= (A.58) 

When the boundary layer trajectory reaches the constraint boundary (i.e. x 2 = 1), a = 0 and the 
expression (A.58) yields a finite final time of 2 units of boundary layer time. Transforming back, 
we obtain an expression for u in terms of x 2 , 


u = V 1 - x 2 

where integration of the original differential equation (A.39) yields. 



(A. 59) 


(A.58) 


V. Conclusions 


In conclusion, state-variable inequality constrained singularly perturbed problems can exhibit 
complex boundary layer phenomenon that are not well understood. The order of the state constraint 
can increase when going from the full order problem to a boundary layer analysis. Because 
discontinuous costate histories can be introduced by the presence of state inequality constraints, a 
direct application of available singular perturbation theory, which requires the state and costate 
histories to be smooth, is not possible. The boundary layer phenomenon associated with such 
problems appear to be finite time. A stability theory for finite time phenomenon, as required to 
construct a suitable approximation for costates appearing in derived feedback control laws, is not 
available at this time. Valentine's transformation can be used to overcome some of these 
difficulties, but at the expense of introducing singular arcs and possibly increased system order. 
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Appendix B 

Permormance Modeling of Hypersonic Vehicles 


I. Introduction 

In order to carry out useful performance studies, trajectory optimization and guidance law 
development for hypersonic transatmospheric vehicles, it is necessary to utilize an accurate 
model of the aerodynamics and propulsion system characteristics. Since actual vehicle design 
is not involved, it is appropriate to utilize simplified models if they can be made to properly 
reflect the actual vehicle performance characteristics. It is also of great benefit to have 
available models that can be used interactively to study the impact of small changes in vehicle 
configuration or propulsion system design. 

In what follows is described a set of simple algorithms devised for use in the present 
research program for the purposes outlined. The computer codes have evolved continuously 
throughout the study. The result is an integrated hypersonic vehicle performance package that 
has many applications beyond those originally envisioned. 

n. Hypersonic Aerodynamic Performance Modeling 

Simple hypersonic aerodynamic theory enables construction of practical and highly accu- 
rate representations of the performance characteristics of realistic hypersonic flight vehicles. 
In this section we review the basic theoretical approach and the implementation of this theory 
in the form of interactive computer software. The basic approach was to make the application 
of the model to a particular airframe conceptual design as simple as possible. Because of the 
interactive nature of the algorithms used, effects of even minor design modifications can be 
immediately assessed in terms of sensitivity parameters such as L/D ratio, overall vehicle drag 
coefficient, and trim moments. 

The models developed have applications that range considerably beyond the ones ad- 
dressed in this report. For example, they are a sufficiently accurate representation of the 
vehicle performance to allow assessment of off-design flight conditions as well as approximate 
stability and control studies. Although the emphasis in the following discussion is on the high- 
Mach number performance modeling, the computer program under development is being set 
up to cover the entire Mach number range from low subsonic to hypersonic speeds. The low 
speed aerodynamic performance models used are not as accurate as those in the hypersonic 
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range, but are sufficiently precise for use in simple performance modeling in which low speed 
flight affects such as the landing or takeoff flight phases are to be included. 

In what follows is given a brief description of the operation of the computer program and 
the basic theory on which it is based. The methods used in applying the theory to a given 
vehicle configuration is reviewed for the benefit of readers unfamiliar with hypersonic aerody- 
namic modeling. 

User Interface 

The computer algorithms were designed to make their application to a given vehicle con- 
figuration as simple as possible. At present, limited access to actual flight vehicle configura- 
tions makes it necessary to work from rather sparse data sets. For example, the vehicle to be 
studied may be defined only by a simple three-view drawing. We have deliberately set out to 
make it possible to work effectively from such data. The configuration is entered into the 
program in a variety of ways. The simplest method allows input in the form of outlines of the 
wing planform, fuselage elevation and planform, body cross-section shapes, and tail surface 
configuration in the form of discrete points. It is not necessary that a large number of outline 
points be used. For example fuselage outline data can consist of as few as ten points in 
elevation and planform with acceptable accuracy. The program allows for variation in 
fuselage cross-section station by station along the axis of symmetry. It also allows for comers 
in the cross-section as often chosen in hypersonic vehicle layouts. 

The configuration data can also be entered by scanning a three-view of the vehicle. Scaling 
is accomplished by selecting points at the nose and tail of the planform. The user then selects 
points interactively by means of a mouse or graphics table. Modifications in geometry can be 
directly implemented in the input process by altering position of control points. Wing and tail 
incidence, fin cant, control surface deflection , and other required information are entered in an 
interactive tabular input window. If insufficient data has been entered to properly define the 
complete configuration, the program warns the operator and indicates what additional informa- 
tion must be specified. 

Interactive Program Mode 

The computer program has been designed to take full advantage of modem computer 
graphical interface technology. Once a vehicle configuration has been implemented as 
described earlier, its attributes can be saved and modified later. At the discretion of the user, 
one or more attributes of the vehicle aerodynamic performance characteristics can be displayed 
simultaneously with the configuration input panel. For example, the lift coefficient vs angle of 
attack, lift/drag ratio, pitching moment vs lift coefficient or other information can be viewed at 
the same time changes in vehicle geometry such as wing area, incidence angle, airfoil shape, or 
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body cross- section shape are being input. This makes it very easy to determine the impact of 
design changes in a direct and graphically useful manner. When the program is linked to one 
of the postprocessing packages such as the one describing trajectory optimization discussed 
elsewhere in the report, it is possible to view effects of configuration modifications directly in 
terms of their impact on selected performance parameters. 

Local Surface Inclination Theory with Blast Wave Corrections 

Several decades of experience have shown that the simplest form of hypersonic flow field 
modeling yields a practical and accurate means for estimating vehicle performance. The 
Newtonian flow model gives an excellent representation for the pressure changes on the 
vehicle surfaces directly in terms of the inclination of the local surface to the ffeestream flow. 
Various modifications can also be applied to correct the pressure distribution for effects of 
strong shock formations at the leading edges of lifting surfaces and tail surfaces and on the 
nose of thefuselage. The blast wave theory is employed for this purpose. 

Simple Newtonian impact theory shows that the pressure coefficient at any point on the 
windward vehicle surface is given by 

C p = 2sin 2 0 (B>1) 

Thus all that is necessary to apply it to a three dimensional vehicle is to set up an algorithm that 
utilizes geometry information to determine whether a given element of the surface is on the 
windward side and to calculate the angle between the freestream velocity vector and the 
surface element. Samples of the computational method used in this program are discussed 
briefly in the following subsections. 

Hypersonic Thin Airfoil Theory 

In some situations, it is sufficiently accurate to represent hypersonic lifting surfaces as flat 
plates. However, in practical situations the need for adequate low-speed aerodynamic charac- 
teristics and surface structure dictates that a cambered airfoil of reasonable thickness be used. 
For example, the NASA vehicle designs used in our computations typically employed airfoils 
of of between 5 and 10% maximum thickness ratio. Thus it is useful to provide means for 
correcting the force calculations for camber and thickness effects. Figure B.l defines the 
required geometry for a typical wing section. 

On the windward side of the airfoil, the camber/thickness function is conveniently de- 
scribed as a functional relationship 

y, = F(x) (B.2) 
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Figure B.l. 



This can be input in the computer program as a series of points taken from a scan of the airfoil 
or by inputing in tabular form. The program also allows the specification of the airfoil profile 
directly in analytical form. The program logic then determines the required calculation module 
from which to compute the wing characteristics. The program also contains provision for 
accounting for wing twist, although none of the vehicle models studied have employed twist. 
The local value of the pressure coefficient becomes 

C p = 2sin 2 (a + l-£) (B.3) 

where a is the vehicle angle of attack measured between the fuselage reference plane (See 
Figure B.2), i is the incidence angle between the wing chord line and the fuselage reference 
plane, and F is the airfoil envelope shape function as described above. This information is then 
used as the basis for determining the normal force coefficient for the airfoil. The result is 
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Figure B.2. Definition of Fuselage Reference Plane 


and from this we find the section lift and drag coefficients. 

(B.5) 
(B.6) 

These coefficients are then used with the wing planform information giving the local values of 
chord length and incidence (for a twisted wing) to determine the force on the local wing 
section. The results for the entire wing are then accumulated. The program has provision for 
displaying the total lift and drag coefficients and the center of pressure location for the three- 
dimensional wing. It also determines contributions to the pitch, yaw, and roll moment 
coefficients. 

Aileron, elevon or flaperon deflection effects are also computed. The user must input the 
desired control surface deflections. Differential elevon deflections are allowed. The program 
senses when the critical surface deflection angle (at which the freestream flow no longer 
impinges on the deflected surface) has been exceeded and properly adjusts the force system. 


2 cos a 
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Hypersonic Lifting Body Theory 


The integrated fuselage/propulsion system provides most of the lift of a typical hypersonic 
aircraft. Since the shapes may be somewhat complex, it is necessary to provide for an adequate 
geometrical representation. Figure B.2 shows the coordinate system used and defines the 
vehicle reference plane. This plane typically coincides with the body centerline as seen in the 
elevation drawing, but the program allows for arbitrary specification of this plane. For 
convenience, the origin of the coordinate system is located at the vehicle nose. The x-axis lies 
along the reference plane in the (assumed) normal plane of symmetry. The y-axis points to the 
left and the z-axis points downward. 

Figure B.3 shows how the various profile curves defining the body shape are represented in the 
program. These curves may be determined by curve fitting of three-view drawings or may be 
input into the program as a table of points. It is not necessary to utilize a large number of 
points. Ten points per profile curve usually provide adequate accuracy unless the body shapes 
are exceptionally complex. 

As shown in Figure B.3, the body cross-section profiles are not required to be continuous 
curves. Comers are allowed as represented by the break points shown in the drawings. The 
fuselage shape is specified in functional form as 

z = f(y) Cross Section Shape 

1 y = g(x) Planform Shape (B.7) 

z = h(x) Fuselage Elevation 


These functions are determined in the program from the input coordinate points and are used to 
compute the unit vector normal to a point on the windward surface. The result is 

n = sin^ji + sin^j + ^os^cos^k (B.8) 


where 


<}>l = tan 



<j >2 = tan 


"df 

,dy 


\ 

J 


(B.9) 

(B.10) 


An area element of the body surface at the same location can be written as 

A 


dS =- 


cos<{>i COS<j>2 


dy 


(B.l 1) 
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Fuselage Elevation 




Figure B.3. Definition of Body Shap e Functions 
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The projected width of the element in the reference plane is A as shown in Figure B.3. The 
inclination angle between the freestream velocity vector and the surface element is given by 


Y = cos 


-ifUco-n 


K 

2 


where the velocity vector is specified in body coordinates as 


(B. 12) 


Uoo = Uoo[(cosP cosa)i-(sinP)j — (cosP sina)k] (B. 13) 


The local pressure coefficient is 


C p = 2 sin 2 [v] 


(B.14) 


The dimensionless normal force vector on the surface element is, in vector form 


dF = -2 sin 2 [\|/] dS n (B.15) 

Integration of this function across the width of the body gives the force vector on an axial 
element of the vehicle. The program then sums all contributions axially to determine the lift, 
drag, sideforce, and moment coefficients for the complete assembly. 

The more detailed representation of the body geometry used in the present version of the 
program greatly improves the agreement between the predicted and measured aerodynamic 
performance. Earlier problems with the vehicle lift curve slope at zero angle of attack have 
been eliminated because the effects of body curvature and surface orientation are now properly 
accommodated in the calculations. 


HI. Scramjet Propulsion Model 

Although the supersonic combustion ramjet concept has been known for over two decades, 
the lack of appropriate unclassified experimental data, cycle analyses, and combustion analy- 
ses requires that we use a simple conceptual model for the purposes of vehicle trajectory opti- 
mization. In what follows is a brief description of this model and the philosophy behind it. 

Conceptually, the SCRAMJET is as simple an airbreathing combustion device as one 
could imagine. In the case of the new family of flight vehicles to use this propulsion concept, 
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the entire underside of the vehicle plays a role in the operation of the system. Figure 1 shows 
the basic configuration. Mechanically, the device can be thought of in terms of three elements. 
These are 

1. Diffuser 

2. Combustor 

3. Expansion nozzle 

Hypersonic vehicle designers attempt to utilize the forward fuselage, strakes, and wings to 
provide the majority of the diffusion. The lower part of the three-dimensional oblique shock 
formed at the leading edges is tailored to the shape of the combustor inlet so that air enters at 
approximately a Mach number of the order of 3 depending on the flight speed. Combustion of 
hydrogen fuel takes place in the duct at supersonic speeds in order to minimize energy losses 
due to dissociation, which would be enormous if the more conventional subsonic ramjet cycle 
were to be used in high speed flight. Liquid hydrogen is the fuel of choice not only because of 
its high energy content, but because it can be made to bum in a supersonic flow due to its wide 
flammability limits and high flame speed. Finally, the combustion products are expanded 
through a nozzle, which, like the diffuser is designed into the contour of the lower fuselage. 

The propulsion system is mostly diffuser and nozzle. While these elements are fairly easy 
to model from the thermodynamic cycle point of view, the aerodynamics are quite complex, 
giving rise to a challenging design problem. Computational fluid dynamics (CFD) numerical 
techniques are being relied upon in conjunction with a new family of hypersonic test facilities 
to yield practical design solutions. Unfortunately, information on the current research is 
classified, so that realistic design data is not available for projects of the type reported here. 

The computational model used here to represent the SCRAMJET propulsion system was 
deliberately designed to be readily updated as new information becomes available. It directly 
accesses a standard atmosphere model (also easily adjustable to provide non-standard operat- 
ing conditions), which simplifies its incorporation into a trajectory optimization program. The 
diffuser and nozzle performance is determined either with standard thermodynamic models or 
by means of optimal design curve fits such as those proposed by Billig. Since information 
concerning recent progress in supersonic combustion was not available, a simple combustor 
model was incorporated. This is a straightforward Rayleigh line calculation. An iterative 
scheme is used to determine the nozzle entrance Mach number, by maintaining the mixture 
ratio at or below the stoichiometric value. No detailed combustion calculations with multi- 
species gases is attempted in the present version of the model although these could be readily 
incorporated as a more definitive model of practical SCRAMJET combustion comes into 
focus. 
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The propulsive drag estimate of Billig was incorporated to account in a simple way for 
some of the frictional losses. No attempt was made to incorporate vehicle integration effects in 
an interactive fashion. Experience with the aerodynamic simulation shows that very small 
vehicle attitude changes take place during equilibrium flight. Therefore in the present state of 
development, no vehicle attitude dependence has been included in the propulsion model. The 
flexibility of the algorithm will make such additions quite easy to incorporate as the need for 
them is established. 
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